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1.  Introduction 


Although  the  analyses  perfonned  by  the  Advanced  Lethality  and  Protection  Analysis  Branch 
(ALPAB)  of  the  Weapons  and  Materials  Research  Directorate,  US  Anny  Research  Laboratory, 
are  extremely  varied  and  diverse,  some  scenarios  arise  with  a  not  unexpected  consistency.  It  is 
not  uncommon  for  ALPAB  analyses  to  require  a  grid  scheme  or  uniform  random  point  picking 
over  a  given  domain.  A  few  general  examples  of  this  are  random  personnel  or  vehicle  targets 
spread  across  an  area  target,  evenly  distributed  sample  points  structured  to  allow  for  unbiased 
spatial  data  collection,  and  space-filling  vehicle  or  projectile  fonnations  (e.g.,  potential  swarming 
applications). 

Point  picking  and  evenly  spaced  point  distributing  over  a  parallelogram  or  parallelepiped  domain 
are  straightforward  and  intuitive.  The  most  basic  example  of  evenly  spaced  point  distributing 
over  a  parallelogram  or  parallelepiped  is  the  set  of  vertices  of  a  Cartesian  grid.  Most  scientists, 
engineers,  and  laypersons  alike  are  familiar  and  comfortable  with  the  concept  of  using  a  Cartesian 
grid  to  evenly  partition  Euclidean  space.  However,  the  process  of  point  picking  or  distributing  over 
other  domains  is  not  always  as  straightforward.  More  specifically,  point  picking  and  especially 
evenly  spaced  point  distributing  over  circular  and  spherical  domains  turns  out  to  be  a  much  more 
difficult  problem,  and  unfortunately,  such  domains  are  regularly  required  or  desired. 

Point  picking  on  the  disc  has  most  recently  and  regularly  been  used  by  ALPAB  analysts  to 
calculate  the  positions  of  individual  personnel  and  vehicle  targets  over  a  circular  area  battlefield. 
One  specific  example  of  such  a  target  was  used  in  the  Modular-Munition,  Dual- Warhead 
Concept  Analysis.1  The  battlefield  target  in  this  analysis  was  supposed  to  loosely  represent  how 
an  assembly  or  staging  consisting  of  both  personnel  and  vehicle  targets  might  look.  This  circular 
area  target  was  also  used  in  subsequent  analyses  such  as  the  Modular  Lethality  Multi-Target 
Defeat  analysis.2  Evenly  spaced  point  distribution  has  most  recently  become  a  tool  in  the 
Modular  Lethality  effort  as  well  and  has  been  applied  to  the  positioning  of  a  set  of  individual 
projectiles  that  are  assumed  to  be  able  to  work  in  a  specific  formation  to  defeat  the  given  target.3 
Evenly  spaced  point  distribution  over  a  disc  or  even  on  or  within  a  sphere  will  most  likely  play  a 
role  in  many  “swarming”  type  concept  analyses  in  the  future. 

The  ALPAB  project  that  finally  spurred  the  collection  of  point  picking  and  distributing 
algorithms  and  the  writing  of  this  report  was  Adaptive  Protection:  the  Dynamic  Threat  Map 
(DTM).  The  objective  of  the  DTM  was  to  develop  a  set  of  algorithms  for  predicting  the  potential 
for  threat  (i.e.,  future  threats)  around  a  vehicle’s  current  location  based  off  of  both  historical  and 
real  time  situational  data.  Part  of  the  algorithm  development  necessitated  the  fast  computation  of 
an  arbitrary  number  of  evenly  spaced  data  collection  rays  emanating  from  the  center  of  mass  of 
the  vehicle.  This  is  exactly  the  problem  of  even  point  distribution  on  the  surface  of  the  sphere 
(hemisphere  in  this  particular  case). 
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These  are  just  a  few  examples  of  the  ALPAB  projects  over  the  past  2  years  that  have 
implemented  algorithms  for  picking  and  distributing  points  on  a  disc  or  sphere.  In  actuality,  the 
author  estimates  that  she  is  asked  about  point  picking  and  distributing  algorithms  at  least  once 
every  other  month.  It  is  with  this  in  mind  that  in  this  report  we  present  a  collection  of  algorithms 
for  picking  and  distributing  points  on  a  disc  and  sphere.  In  computer  simulation,  the  need  to 
arrange  points  on  a  disc  or  on  or  within  a  sphere  commonly  falls  into  2  categories:  1)  random 
point  picking  and  2)  evenly  spaced  point  distributing. 

Random  point  picking,  or  just  point  picking,  is  defined  in  this  report  as  the  process  of  picking 
random  points  within  the  region  of  interest  such  that  the  density  of  points  within  arbitrary 
neighborhoods  is  expected  to  be  unifonn  over  the  entire  region  as  the  number  of  points,  n,  goes 
to  infinity.  Random  point  picking  is  very  useful  when  an  unstructured  (i.e.,  random)  point 
distribution  is  desired.  However,  for  small  values  of  n,  random  point  picking  does  not  guarantee 
uniform  spatial  distribution  of  points  over  the  region. 

Evenly  spaced  point  distributing  is  defined  in  this  report  as  the  process  of  distributing  n  points 
over  the  surface  of  the  disc  or  sphere  so  as  to  produce  the  most  evenly  spaced  distribution  of 
points  possible.  To  consistently  produce  unifonn  and  evenly  spaced  distributions  of  points  when 
n  is  small,  evenly  spaced  point  distributing  techniques  often  produce  structured,  grid-like  point 
patterns. 

This  report  does  not  attempt  to  document  an  exhaustive  collection  of  point  picking  and 
distributing  algorithms  but  instead  focuses  on  the  small  handful  that  the  author  has  found  to  be 
most  useful  for  her  purposes.  All  calculations  are  presented  for  a  disc  or  sphere  centered  at  the 
origin. 


2.  Point  Picking  on  (and  around)  the  Disc 

2.1  Elimination  Method 

If  runtime  is  not  an  issue,  an  easy  method  for  point  picking  on  a  disc  of  radius  r  is  the 
elimination  method.  In  the  elimination  method,  the  x-  and  y-coordinates  of  each  point  are  each 
drawn  from  a  uniform  [— r,  r]  distribution,  (i.e.,  from  the  smallest  square  containing  the  disc). 
Each  point  is  then  either  accepted  or  rejected  depending  on  whether  or  not  it  lies  within  the  disc 
(Fig.  1).  This  method  is  simple  to  implement,  but  as  the  disc  only  occupies  approximately  79% 
of  the  area  within  this  square  domain,  approximately  21%  of  the  points  picked  will  be  rejected 
(Fig.  1).  For  this  reason,  it  is  often  desired  to  implement  a  direct  method  that  does  not  rely  on 
picking  and  eliminating  points. 
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Point  Picking  on  the  Disc:  Elimination  Method 


^—Minimal  square 
- Disc  of  radius  r 


Accepted  points 


Fig.  1  Elimination  method  for  point  picking  on  the  disc.  Points  are  generated  over  the 
smallest  square  region  containing  the  disc  and  then  are  rejected  if  they  do  not  lie 
within  the  disc.  Only  about  79%  of  points  will  be  accepted  with  approximately 
21%  being  rejected. 


2.2  Direct  Method:  The  Most  Common  Mistake 

When  point  picking  over  a  rectangular  domain  [a,  b]  x  [ c ,  d]  c  R2,  most  people  will  quickly 
realize  that  the  most  straightforward  technique  is  to  draw  an  x-coordinate  from  a  uniform  [a,  b] 
distribution,  and  a  y-coordinate  from  a  unifonn  [ c ,  d]  distribution.  The  resulting  points  will 
exhibit  a  unifonn  spatial  distribution  over  the  domain  (Fig.  2). 
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Point  Picking  on  a  Rectangle 


Fig.  2  Point  picking  on  a  rectangle;  x-  and  y-coordinates  are  drawn  independently  from  uniform 
distributions  over  their  respective  domains  and  yield  a  uniform  spatial  distribution  over 
the  rectangle 


The  logical  extension  of  this  technique  to  point  picking  over  a  circular  domain 

S 1  ■—  {(x,y)  £l2  |  x2  +  y2  <  r2}  (1) 

is  to  switch  to  polar  coordinates  over  the  domain  [p,  0]  Q  [0,  r]  x  [0,27r] .  We  can  then  draw  p 
from  a  uniform  [0,  r]  distribution,  and  0  from  a  uniform  [0,27r]  distribution,  and  use  the  standard 
polar  to  Cartesian  transformations 


x  —  p  cos  6 

(2) 

y  —  p  sin  6 

(3) 

to  get  back  to  the  desired  Cartesian  coordinates.  Although  the  resulting  points  will  have  p-  and 
6 -coordinates  that  are  independently  uniformly  distributed,  it  quickly  becomes  apparently  that 
the  points  themselves  are  “clumped”  near  the  origin  and  far  from  uniformly  distributed  over  the 
disc  (Fig.  3). 
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Point  Picking  on  the  Disc 
Direct  Method  -  Incorrect 


Fig.  3  Point  picking  on  the  disc;  drawing  both  p-  and  61- 

coordinates  from  independent,  uniform  distributions 
over  their  respective  domains  yields  a  point  distribution 
that  is  clumped  around  the  center  of  the  disc 

There  is  a  simple  explanation  for  this  outcome.  First,  recall  that  the  area  of  the  sector  subtended 
by  a  central  angle,  <p,  of  a  circle  of  radius  r  (Fig.  4)  is  given  by  the  equation 

A  =  ^r2  <p.  (4) 

Therefore,  the  area  of  a  polar  rectangle  about  a  point,/?,  as  shown  in  Fig.  5  can  be  expressed  as 

dA  -  A2  -  Ax  =  \pl  dd  -  ipf  dd  =  i  (p|  -  pi)  dd.  (5) 

This  can  be  expanded  to 

dA  =  \  (p2  +  Pi)(p2  -  Pi)  dd  =  i  (p2  +  px)  dp  dd  =  p  dp  dd,  (6) 

where 

P  =  ^(p2+Pi)-  (V) 

Equation  6  is  commonly  referred  to  as  the  area  element  and  should  be  familiar  to  anyone  who 
has  evaluated  a  double  integral  in  polar  coordinates. 

So  when  point  picking  on  the  disc,  although  the  density  of  p-coordinates  over  any  change  in  p 
(i.e.,  over  any  dp)  and  the  density  of  0-coordinates  over  any  change  in  0  (i.e.,  over  any  dd) 
remain  the  same,  the  area  of  the  resulting  polar  rectangle  (i.e.,  the  area  element)  is  given  by 

dA  —  p  dp  dd.  (8) 
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Area  of  a  Sector 


Fig.  4  Area  of  a  sector  of  a  circle 


Polar  Rectangle 


Fig.  5  Example  of  a  polar  rectangle 
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This  means  that  for  any  2  polar  rectangles  with  equal  dps  and  dOs,  the  number  of  points  within 
the  polar  rectangles  will  be  equal;  however,  the  area,  and  therefore  density  of  points  within  each 
polar  rectangle  will  directly  depend  on  the  p-coordinate  of  that  particular  polar  rectangle  (Fig.  6). 


Polar  Rectangles  with  Unequal  Point  Densities 


Fig.  6  Two  area  elements  with  equal  radial  lengths  (dp)  and 
angular  components  (dd)  but  different  areas 


2.3  Direct  Method:  A  Simple  Transform 

Luckily,  there  is  a  rather  intuitive  solution  to  the  “clumping  problem”.  We  need  only  to  find  the 
appropriate  coordinate  transform  to  remove  the  radial  coordinate  dependence  from  the  area  of 
the  polar  rectangle.  To  be  exact,  we  are  looking  for  a  coordinate  transfonn  of  the  form 

u  =  f(p ')»  f-  [0,r]  ->  [a,b],  f^C1,  (9) 

such  that 


du  =  p  dp. 

Then  the  area  element  in  Eq.  8  would  become 

dA  —  p  dp  d6  —  du  dO. 

Integrating  Eq.  1 1  yields  the  desired  transformation 

u  =  ^p2,  u\  [0,r]  ->  0,ir2]. 
For  convenience,  we  multiply  by  2  and  adjust  the  interval  to  yield 


(10) 

(11) 

(12) 
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u  =  p2,  it:  [0,r]  ->  [0 ,r2]. 


(13) 


Finally,  to  make  use  of  the  transfonn,  we  solve  for  p  to  yield 

p  =  yfu,  p\  [0,r2]  -»  [0,r]. 


(14) 


With  this  coordinate  transform,  the  direct  method  for  point  picking  over  a  disc  of  radius  r 
reduces  to  drawing  it  from  a  unifonn  [0,r2]  distribution,  6  from  aunifonn  [0, 2ttJ  distribution, 
and  applying  the  transformations 


and 


x  =  Vwcos  6 

(15) 

y  =  yfu  sin  6. 

(16) 

The  resulting  points  no  longer  exhibit  the  clumped  distribution  near  the  origin  (Fig.  7),  and  it  can 
be  easily  verified  that  the  points  are  now  unifonnly  distributed  across  the  disc  (left  to  the  reader). 


Direct  Point  Picking  on  the  Disc 


Incorrect 


Correct 


Fig.  7  Side  by  side  comparison  of  the  most  common  mistake  and  correct  method  for  directly 
picking  points  on  the  disc 

2.4  Points  on  the  Circle 

One  application,  albeit  trivial,  of  the  direct  point  picking  method  presented  in  Section  2.3  is 
picking  uniformly  distributed  points  along  the  boundary  of  the  disc,  (i.e.,  around  a  circle).  The 
boundary  of  the  disc  is  the  set  of  points 

dS1  {(x,y)  G  R2  |  x2  +  y2  —  r2},  (17) 


which  is  equivalent  to  the  set 
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H1  :=  { 6  G  [0,27t]  |  x  —  r  cos  6,  y  —  r  sin  6}. 


(18) 


Therefore,  the  direct  method  for  point  picking  around  a  circle  of  r  reduces  to  drawing  9  from  a 
uniform  [0,27r]  distribution  and  applying  the  transfonnations 


and 


x  —  r  cos  6 


(19) 


y  —  r  sin  6 . 


(20) 


3.  Point  Picking  in  (and  on)  the  Sphere 


3.1  Elimination  Method 

If  runtime  is  not  an  issue,  an  easy  method  for  point  picking  within  a  sphere  of  radius  r  is  the 
elimination  method.  In  the  elimination  method,  the  x-,  y-,  and  z-coordinates  of  each  point  are 
each  drawn  from  a  uniform  [— r,  r ]  distribution,  i.e.,  from  the  smallest  cube  containing  the 
sphere.  Each  point  is  then  either  accepted  or  rejected  depending  on  whether  or  not  it  lies  within 
the  sphere  (Fig.  8).  This  method  is  simple  to  implement,  but  as  the  sphere  only  occupies 
approximately  53%  of  the  volume  within  this  cube  domain,  approximately  47%  of  the  points 
picked  will  be  rejected  (Fig.  8).  For  this  reason,  it  is  often  desired  to  implement  a  direct  method 
that  does  not  rely  on  picking  and  eliminating  points. 

Point  Picking  in  the  Sphere:  Elimitation  Method 


53%  inside 
the  sphere 


47%  outside  the  sphere 


Fig.  8  Elimination  method  for  point  picking  within  the  sphere.  Points  are  generated  within 
the  smallest  cubic  region  containing  the  sphere  and  then  are  rejected  if  they  do  not  lie 
within  the  sphere.  Only  about  53%  of  points  will  be  accepted  with  approximately  47% 
being  rejected.  For  visual  clarity,  only  points  in  the  “back  half’  of  the  region  are 
shown  in  the  graph  on  the  left. 
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3.2  Direct  Method:  The  Most  Common  Mistake 


As  in  the  2-dimensional  (2-D)  analogue,  when  point  picking  over  the  domain 

S2  {(x,y,z)  £R3  |  x2  +  y2  +  z2  <  r2},  (21) 

most  people  will  switch  to  spherical  coordinates  over  the  domain  [p,  6,  cp]  Q  [0,  r]  x  [0,27rJ  x  [0,  tt], 
attempt  to  draw  each  coordinate  from  a  unifonn  distribution  over  its  respective  interval,  and  transfonn 
back  to  Cartesian  coordinates  using  the  standard  transformations 


x  —  p  cos  6  sirup 

(22) 

y  —  p  sind  sin<p 

(23) 

z  —  p  cos  cp. 

(24) 

Unfortunately,  this  once  again  produces  clumped  points;  this  time  they  appear  to  be  clumped 
along  the  z-axis  and  most  tightly  at  the  origin  (Fig.  9). 

Incorrect  Point  Picking  in  the  Sphere 


Side  View  Top  View 


Fig.  9  Point  picking  within  the  sphere;  drawing  p-,9-,  and  (^-coordinates  from  independent,  uniform 
distributions  over  their  respective  domains  yields  a  point  distribution  that  is  clumped  along  the 
z-axis  and  most  tightly  near  the  center  of  the  sphere. 


Again,  although  the  density  of  p-coordi  nates  over  any  change  in  p  (i.e.,  over  any  dp),  and  the 
density  of  0-coordinates  over  any  change  in  6  (i.e.,  over  any  dd ),  and  the  density  of  <p- 
coordinates  over  any  change  in  cp  (i.e.,  over  any  d(p)  remain  the  same,  the  volume  of  the 
resulting  spherical  quadrilateral  (Fig.  10)  is  given  by 

dV  —  p2  sin  (p  dp  d6  dcp.  (25) 
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Equation  25  is  commonly  referred  to  as  the  volume  element  and  should  be  familiar  to  anyone 
who  has  evaluated  a  triple  integral  in  spherical  coordinates. 


Spherical  Quadrilateral 


Fig.  10  Example  of  a  spherical  quadrilateral 


This  means  that,  for  any  2  spherical  quadrilaterals  with  equal  dp s,  dcp s,  and  dd s,  the  number  of 
points  within  the  spherical  quadrilaterals  will  be  the  same.  However,  the  volume,  and  therefore 
density,  of  points  within  each  spherical  quadrilateral  will  depend  on  both  the  p-  and  (p- 
coordinates  of  the  spherical  quadrilateral. 

3.3  Direct  Method:  A  Simple  Transform 

Luckily,  the  solution  to  the  clumping  problem  in  3-D  is  similar  to  that  in  2-D.  We  need  only  to 
find  the  appropriate  coordinate  transforms  to  remove  the  p-  and  (^-coordinate  dependences  from 
the  volume  of  the  spherical  quadrilateral.  In  this  case,  we  are  looking  for  a  coordinate  transfonns 
of  the  form 


u  —  f  (p),  /:  [0,r]  ->  [a,b],  f  E  C1 
v  =  g{<p ),  g\ [0,7r]  ->  [c,d],  g  6  C1 


(26) 


(27) 


such  that 


du  =  p2  dp 


(28) 


and 


dv  —  sirup  dq). 


(29) 
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Then  the  volume  element  in  Eq.  25  would  become 


dV  =  p 2  sin  q)  dp  dcp  d9  =  du  dv  dO.  (30) 

Note  the  change  in  order  of  integration.  This  is  possible  as  we  will  now  be  integrating  over  a 
constant  function  with  pairwise  independent  limits  of  integration. 

Integrating  Eqs.  28  and  29  yields  the  desired  transfonns 


1  3 

U=  3  P  ’ 


u:[0,r]-»  0,^r3 


(31) 


and 


v  =  —co scp,  v:  [0, 7r]  -*  [— 1,1]  .  (32) 

For  convenience,  we  can  drop  the  constant  coefficients  and  adjust  the  intervals  to  yield 

u  =  p3,  it:  [0,  r]  [0,  r3]  (33) 

and 

v—  cos  (p,  v:  [0, 7r]  -*  [— 1,1].  (34) 

Finally,  in  order  to  make  use  of  these  transfonns,  we  take  the  inverse  of  Eq.  33  to  yield 

p=Vu,  P-  [0,r3]  -»  [0,r],  (35) 

and  make  the  following  observation  regarding  Eq.  34: 

sin<p  —  -J 1  —  cos2  cp  —  Vl  —  v2,  v  G  [—1,1].  (36) 

With  these  coordinate  transfonns,  the  direct  method  for  point  picking  within  a  sphere  reduces  to 
drawing  it  from  a  uniform  [0,  r3J  distribution,  v  from  a  unifonn  [—1, 1]  distribution  and  6  from 
a  unifonn  [0,27r]  distribution,  and  applying  the  transformations 

x  =  Vu  Vl  —  v2  cos  9  ,  (37) 

y  —  \fu  Vl  —  v2  sin  6  ,  (38) 

and 


z  =  Mu  v. 


(39) 


The  resulting  points  no  longer  exhibit  the  clumped  distribution  near  the  origin  (Figs.  1 1  and  12), 
and  it  can  be  easily  verified  that  the  points  are  now  uniformly  distributed  in  x,  y,  and  z  (left  to  the 
reader). 
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Correct  Point  Picking  in  the  Sphere 


Side  View 


Top  View 


Fig.  11  Point  distribution  using  the  correct  parameter  domains  and  transformations.  Notice  the  absence 
of  the  clumping  problem. 


Direct  Point  Picking  in  the  Sphere 


Incorrect  Correct 


Fig.  12  Side  by  side  comparison  of  the  most  common  mistake  and  correct  method  for  directly  picking 
points  within  the  sphere 
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3.4  Points  on  the  Surface  of  the  Sphere 

Less  trivial  than  its  2-D  counterpart,  one  extension  of  the  direct  point  picking  method  presented 
in  Section  3.3  is  picking  uniformly  distributed  points  on  the  surface  of  the  sphere  (Fig.  13). 
There  are  2  common  methods  for  achieving  this. 

Point  Picking  on  the  Sphere 


Fig.  13  Point  picking  on  the  sphere.  Notice  that  the  points  appear  to  cluster  around  the  image  boundary 
in  both  views.  This  is  due  to  the  2-D  visualization  of  the  points  on  the  surface  of  the  3-D  sphere. 


The  first  is  the  projection  method.  In  the  projection  method,  the  desired  n  points  are  first 
generated  within  the  sphere  using  any  method  desired.  The  points  are  then  projected  onto  the 
surface  of  the  sphere  by  normalizing  their  Cartesian  coordinates.  This  method  is  straightforward, 
intuitive,  and  simple  to  implement,  but  involves  calculating  the  Euclidean  norm  of  each  point 
and  dividing  the  x-,  y-,  and  z-coordinates  of  each  point  by  its  respective  norm. 


The  second  method  does  not  require  any  nonn  calculations  or  divisions  and  instead  makes  use  of 
the  equations  from  Section  3.3.  The  user  need  only  note  that  the  surface  of  the  sphere  is  the  set  of 
points 


dS 2  {(x,  y,  z)  £  R3  |  x2  +  y2  +  z2  =  r2}, 

which  is  equivalent  to  the  set 

x  =  r  sin<p  cos  0' 
>2  •—  1  fa  e  rn  9-n-t  v  rn  ttI  y  =  r  cos  cp  cos  6  [. 

z  —  r  cos  (p 


ft2  6  [0,27t]  x  [0, 7r] 


(40) 


(41) 
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Therefore,  the  direct  method  for  point  picking  on  the  surface  of  a  sphere  of  r  reduces  to  drawing 
v  from  a  uniform  [—1, 1]  distribution  and  9  from  a  unifonn  [0,27t]  distribution  and  applying  the 
transformations 


and 


x  =  rVl  —  v2  cos  9  , 
y  =  rVl  —  v2  sin  9  , 


z  -  rv. 


(42) 

(43) 

(44) 


4.  Point  Picking  on  an  )V-Sphere 


An  alternate  method  for  point  picking  on  the  surface  of  the  sphere  is  presented  here  as  a  general 
method  that  extends  to  higher  dimensions.  This  general  method,  which  we  will  call  Muller’s 
method,  is  presented  without  proof,  and  the  reader  is  directed  to  the  specified  references  for 
insight  and  proof  for  the  validity  of  this  method. 

Muller’s  method4,5  is  a  2-step  process  for  uniformly  generating  points  on  the  boundary  of  the  N- 
sphere, 

SN  :=  (x  =  (xltx2t -,XN+1)  £  Rw+1  |  Ek=i  Xfc  =  r2}.  (45) 

Step  1:  Generate  N independent  standard  nonnal  random  variables,  yfe,  k  =  1, ... ,  N  +  1. 

Step  2:  Construct  the  point 

x  =  (xllx2,  ...,xN+1)  (46) 

where 

Xk  =  (r/\\y\\)yk,  k  =  1 . N  +  1  (47) 

and 

llyll  =  (iKlyk)1'2-  (48) 

Figure  14  shows  points  generated  by  Muller’s  method  for  N  —  2.  It  can  be  easily  verified  that 
these  points  are  uniformly  distributed  over  the  surface  of  the  sphere  (left  to  the  reader). 
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Muller's  Method  for  Point  Picking  on  the  Sphere 


Side  View 


Top  View 


Fig.  14  Point  picking  on  the  sphere  using  Muller’s  Method.  Notice  that  the  points  appear  to  cluster 

around  the  image  boundary  in  both  views.  This  is  due  to  the  2-D  visualization  of  the  points  on 
the  surface  of  the  3-D  sphere.  This  distribution  appears  to  be  equivalent  to  that  generated  with 
the  direct  method  presented  in  Section  3.4.  The  distributions  are  truly  equivalent,  and  the  proof 
is  left  to  the  reader. 


5.  Evenly  Spaced  Points  on  a  Disc 


Algorithms  presented  in  this  section  for  evenly  distributing  points  across  a  disc  are  not  optimal 
solutions.  In  fact,  the  concept  of  what  it  means  to  be  an  evenly  spaced  distribution  is  not  going  to 
be  formally  defined.  Instead,  we  offer  a  “visual  definition”  in  the  fonn  of  2  graphs  in  Fig.  15. 

The  points  in  the  graph  on  the  left  are  not  evenly  spaced,  while  those  in  the  graph  on  the  right  are 
evenly  spaced. 


16 


Visual  Definition  of  Even  Spacing 


Not  Evenly  Spaced 


Evenly  Spaced 


Fig.  15  Visual  definition  of  what  it  means  for  points  to  be  “evenly  spaced”.  The  points  within  the  circle 
on  the  left  are  not  evenly  spaced.  The  points  within  the  circle  on  the  right  are  evenly  spaced. 

There  are  many  methods  that  will  produce  more  evenly  spaced  points  than  those  presented  in  this 
report,  but  those  methods  often  come  at  a  significant  computational  cost.  Many  of  these  methods 
involve  treating  a  set  of  points  as  particles  confined  to  the  disc,  and  minimizing  the  forces 
between  the  particles.  Especially  as  the  number  of  particles  increases,  such  an  optimization 
problem  can  become  computationally  impractical.  Instead,  we  look  for  approximate  methods  that 
only  involve  algebraic  expressions. 

When  attempting  to  evenly  distribute  n  points  across  the  surface  of  a  disc  of  radius  r  using  only 
algebraic  equations,  a  common  first  idea  is  to  employ  some  form  of  spiral  algorithm.  In  fact, 
there  are  entire  families  of  spirals  that  work  very  well  for  this  purpose.  One  of  the  most  useful 
families  of  spirals  are  the  Archimedean  spirals.  In  polar  coordinates,  Archimedean  spirals  all 
have  the  form 

p  —  a  +  bO 1^c,  a,  b  E  R,  c  G  R\{0}.  (49) 

The  following  sections  describe  how  to  use  2  subfamilies  of  Archimedean  spirals,  along  with  the 
help  of  the  spiral  of  Theodorus,  to  produce  evenly  spaced  point  distributions  on  a  disc  of 
radius  r. 

5.1  Single  Spiral  Method:  Archimedes  and  Theodorus 

The  idea  behind  the  single  spiral  method  is  to  construct  a  sequence  of  n  points  along  a  single 
spiral  that  starts  at  the  origin,  winds  out  toward  the  disc’s  boundary,  and  has  the  following 
properties: 
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1 .  The  distance  between  any  2  successive  points  along  the  spiral  is  constant. 

2.  The  distance  between  any  2  successive  windings  of  the  spiral  is  constant. 

3.  The  distance  between  any  2  successive  points  along  the  spiral  is  as  close  as  possible  to  the 

distance  between  any  2  successive  windings  of  the  spiral. 

Archimedes’s  spiral  (Fig.  16)  is  the  name  for  the  case  when  c  —  1  in  Eq.  49.  Archimedes’s  spiral 
exhibits  the  useful  characteristic  that  the  distance  between  any  2  successive  windings  along  any 
ray  emanating  from  the  origin  is  constant  (Fig.  17).  In  fact,  this  distance  is  known  exactly,  and 
has  the  equation 

d  =  2nb.  (50) 


Achimedes's  Spiral 
Archimedean:  c  =  1 


Fig.  16  Archimedes’s  spiral  is  the  Archimedean  spiral  given  by 
p  =  a  +  b6 ,  where  a,  b  E  E 
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Achimedes's  Spiral 
Archimedean:  c  =  1 


Fig.  17  A  useful  characteristic  of  Archimede’s  spiral  is  that 
the  distance  between  any  2  successive  windings  is 
constant.  Only  the  positive  arm  is  shown;  only  one 
arm  of  the  spiral  was  needed  to  derive  the  single 
spiral  method,  and  the  positive  arm  was  chosen 
arbitrarily. 

This  characteristic  satisfies  property  (2)  and  makes  Archimedes’s  spiral  an  ideal  contender  for 
constructing  our  single  spiral  point  distribution.  The  next  step  is  to  pick  a  and  b,  and  either  p  or 
6  in  such  a  way  that  the  distance  between  any  2  successive  points  along  the  spiral  is  constant. 

The  spiral  of  Theodorus  is  a  spiral  that  was  first  constructed  by  the  Greek  mathematician 
Theodorus  of  Cyrene  in  the  5th  century  BCE.6  The  spiral  of  Theodorus  starts  with  an  isosceles 
right  triangle  and  is  formed  by  a  sequence  of  contiguous  right  triangles  such  that  hypotenuse  of 
each  triangle  become  the  longer  cathetus  (i.e.,  the  longer  nonhypotenuse  side)  of  each  successive 
triangle  (Fig.  18).  A  very  useful  characteristic  of  the  spiral  of  Theodorus  is  that  the  Euclidean 
distance  between  any  2  successive  vertices  of  the  spiral  is  constant. 
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Spiral  of  Theodorus 


Fig.  18  The  spiral  of  Theodorus,  also  known  as  the 
square  root  spiral,  Einstein  spiral,  or 
Pythagorean  spiral 

Similar  to  Archimedes’s  spiral,  the  spiral  of  Theodorus  exhibits  the  characteristic  that  the 
distance  between  any  2  successive  windings  along  any  ray  emanating  from  the  origin  is  nearly 
constant  (Fig.  19).  Although,  this  does  not  exactly  satisfy  property  (2),  the  distance  between  any 
2  successive  winding  monotonically  approaches  the  value  of  n  from  above  very  rapidly. 
Furthermore,  as  n  ->  oo,  the  spiral  of  Theodorus  happens  to  be  an  extremely  close  approximation 
to  Archimedes’s  spiral  with  the  equation 

P  =  \W-  !)•  (51) 
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Spiral  of  Theodorus 


Fig.  19  Similar  to  Archimedes’s  spiral,  the  distance 
between  successive  windings  of  the  spiral  of 
Theodorus  is  nearly  constant 


These  2  spirals  together  provide  all  the  tools  we  need  to  construct  a  single  spiral  of  n  points  that 
satisfies  our  properties  (1),  (2),  and  (3). 

Starting  with  the  sequence  of  hypotenuses  from  the  spiral  of  Theodorus  as  the  radial  coordinates, 
we  fonn  the  sequence  of  points 

(Vk  =  ( Pk>  Qk)  =  (Vfc,  °kj)k=i  (52) 


At  this  point,  the  angular  components  in  the  sequence  are  completely  unknown.  We  then  break 
the  structure  of  the  spiral  of  Theodorus  by  adjusting  the  “winding  constant”,  the  coefficient  of  6 
in  Eq.  49,  until  property  (3)  is  satisfied.  In  other  words,  we  are  going  to  look  for  Archimedes’s 
spiral  such  that  the  distance  between  consecutive  windings  is  as  close  as  possible  to  the  distance 
between  points  with  consecutive  p-values  from  Eq.  52.  More  precisely,  we  are  looking  for  an 
equation  of  the  form 

9  =\p,  b  £  R+  (53) 

such  that 


distance(pk,  pk+1 )  =  2nb. 
Expanding  Eq.  54  in  polar  coordinates  yields 


(54) 


V Pk 2  +  Pk+ 12  2  Pk  pk+i 


cos(0fc+1  -  0fc)  =  2  nb. 


(55) 
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(56) 


Substituting  Eq.  53  into  Eq.  55  and  simplifying  yields 

Pk2  +  Pk+i2  ~2pk  Pk+i  cos  (**p)  =  4  t z2b2. 

Finally,  substituting  Eq.  52  into  Eq.  56  and  simplifying  yields  the  relationship 

2k  +  1  —  2  Vk  yjk  +  1  cos  =  4  n2b2 .  (57) 

Although  it  may  not  be  apparent,  taking  the  limit  as  k  -»  n  and  n  -»  oo  of  the  left  hand  side  of 
the  equation  reveals  the  solution  (see  the  Appendix  for  full  limit  evaluation) 

lim  lim  2k  +  1  —  2  V/r  V /e  +  1  cos  f^k+1  =  -3—  (58) 

n->oo  fc->n  V  ft  /  4b2 

Therefore,  as  k  ->  n  and  n  -»  oo,  Eq.  57  becomes 

-^  =  4  7r2b2,  (59) 

and  we  are  left  with  the  desired  turning  parameter 

b=irn  <60> 

Putting  Eqs.  52,  53,  and  60  together  produce  the  sequence  of  points 

(Pk  =  ( Pk.  9k)  =  (Vfc,  2V7rk))fc=i  ,  (61) 

which  satisfy  all  3  desired  properties.  However,  these  points  cover  a  disc  of  radius  Vn  instead  of 
a  disc  of  the  desired  radius  r.  The  final  step  is  to  introduce  a  radial  scale  factor 

pk  =  tVk ,  t  G  E+  (62) 

and  impose  the  radial  constraint 


Pn  =  r. 

T 

Eqs.  62  and  63  uniquely  detennine  a  value  t  =  j=,  and  produce  the  sequence  of  p-valucs 


(63) 


(64) 


Since  scaling  the  radial  component  of  a  set  of  points  in  polar  coordinates  by  a  constant  value 
scales  the  pair-wise  distances  between  points  by  that  value,  scaling  the  p-values  this  way  not 
only  scales  the  distance  between  successive  windings  by  a  constant  value,  but  it  also  scales  the 
distance  between  successive  points  on  the  spiral  by  the  same  amount.  So  Eqs.  61  and  64  produce 
the  desired  sequence  of  points 
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n 


(65) 


Pk  =  ( Pk >  9k)  =  r 


^  • 
J  '  k=l 


Figure  20  shows  10,000  evenly  spaced  points  generated  by  the  single  spiral  method. 


Evenly  Spaced  Points 
Single  Spiral  Method 


Fig.  20  Ten  thousand  evenly  spaced  points  on  a  disc 
generated  by  the  single  spiral  method 

5.2  Single  Spiral  Method  Revisited:  Fixing  the  Center  Point 

Although  the  sequence  of  points  formed  from  Eq.  65  produce  a  fairly  unifomi  covering  of  the 
disc,  the  reader  may  notice  in  practice  that  this  method  always  leaves  a  bit  of  a  blank  space  near 
the  center  of  the  of  disc.  The  following  equations  are  provided  as  “patch”  to  fill  this  space.  This 
“patch”  was  developed  solely  by  trial  and  error,  but  the  method  appears  to  work  well  in  practice 
for  any  value  of  n. 


Figure  2 1  shows  the  patch  for  2  values  of  n. 
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Single  Spiral  Method  with  "Patched"  Center 


n  =  100 


n  =  1000 


Fig.  21  Single  spiral  method  with  “center”  point  adjusted  to  better  fill  the  space 

From  this  point  on,  any  reference  to  the  single  spiral  method  refers  to  this  “patched”  single  spiral 
method. 

5.3  Vogel’s  Method:  Fermat 

Quite  often,  nature  does  it  best.  Consider  the  head  of  the  sunflower  (Fig.  22),  the  seeds  of  which 
are  evenly  distributed  over  the  “disc”.  Extensive  research  regarding  the  pattern  formed  by  the 
seeds  in  the  head  of  a  sunflower  exists  for  the  interested  reader,  but  for  the  purpose  of  this  report, 
we  acknowledge  the  spiral  phenomena  and  focus  on  an  implementation  of  some  work  done  by 
Helmut  Vogel.7  Vogel  suggested  the  use  of  Fermat’s  spiral,8  the  case  when  c  —  2  in  Eq.  49 
instead  of  the  commonly  used  logarithmic  spiral  to  model  the  pattern  of  the  seeds  on  the  head  of 
a  sunflower  (Fig.  23). 


24 


Fig.  22  Spirals  can  be  found  all  over  the  natural  world  such  as  in  the  arrangement  of  seeds  on  the  head  of  a 
sunflower 


Fermat's  Spiral 
Archimedean:  c  =  2 


Logarithmic  Spiral 


Fig.  23  Both  logarithmic  and  Fermat’s  spirals  can  be  used  to  model  the  observed  spiral  of  seeds  on  the 
head  of  a  sunflower 
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Vogel’s  model  then  has  the  fonn 


and 


6V  =  k  *  8 


Pk  =  tyfk,  G 


(68) 


(69) 


1+V5 


where  8  is  the  golden  angle.  The  golden  angle  is  related  to  the  golden  ratio,  <p  —  — — ,  by  the 
equation 


8  —  2n  ^1  —  —  )  =  7t(3  —  a/5), 


(70) 


where  the  golden  ratio  is  the  limit  of  the  ratios  of  successive  terms  in  the  Fibonacci  sequence. 

Again  imposing  the  radial  constraint  in  Eq.  63,  we  arrive  at  the  final  form  of  what  we  will  call 
Vogel’s  method9 


Pk  =  ( pk.  Ok)  =  rjf .  rc(3  -  V5 )(k-  1) 


(71) 


k=  1 


Figure  24  shows  10,000  evenly  spaced  points  generated  by  Vogel’s  method. 

Evenly  Spaced  Points 
Vogel's  Method 


Fig.  24  Ten  thousand  evenly  spaced  points  on  a  disc  generated 
by  Vogel’s  method 
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5.4  Spiral  Method  Comparison 

The  question  must  now  be  asked:  Is  one  spiral  method  better  than  the  other?  Unfortunately,  there 
is  no  clear  answer,  especially  as  we  have  avoided  formally  defining  what  it  means  to  be  a  “good” 
method  for  generating  evenly  spaced  points.  One  possible  metric  that  can  be  used  to  compare 
these  (and  other)  methods  for  generating  evenly  spaced  points  is  the  standard  deviation  in  the 
size  of  the  cells  of  the  Voronoi  diagram  of  the  points.  A  Voronoi  diagram  is  a  partitioning  of  the 
plane  into  convex  polygons  with  respect  to  a  set  of  vertices  (one  vertex  per  cell)  such  that  any 
point  within  a  given  cell  is  closer  to  that  cell’s  vertex  than  it  is  to  any  other  vertex.  So  while  the 
average  area  of  any  partitioning  of  S1  into  n  regions  will  always  be  n  *  r2/n,  the  standard 
deviation  of  the  area  of  regions  can  provide  insight  into  how  much  the  sizes  of  the  regions  vary. 

Figure  25  shows  Voronoi  diagrams  and  standard  deviation  of  cell  area  for  both  spiral  methods 
for  500  points.  Although  both  methods  produce  distributions  of  points  that  appear  to  be  fairly 
evenly  spaced,  the  standard  deviation  of  the  area  of  Voronoi  cells  of  Vogel’s  method  is  slightly 
more  than  28%  less  than  that  of  the  single  spiral  method.  This  indicates  that  Vogel’s  method 
produces  a  more  evenly  spaced  set  of  points  as  measured  by  the  stated  metric. 

Spiral  Method  Comparison 
n  =  500 


Single  Spiral  Method 

Voronoi  Cell  Standard  Deviation  =  7.313272e-04 
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Vogel's  Method 

Voronoi  Cell  Standard  Deviation  =  5.238220e-04 
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Fig.  25  Voronoi  diagram  comparison  of  single  spiral  and  Vogel’s  method  for  generating  evenly  spaced 
points  on  a  disc 
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6.  Evenly  Spaced  Point  Distribution  on  a  Sphere 


The  problem  of  “evenly  spacing”  n  points  on  the  surface  of  the  sphere  was  first  proposed  by 
JJ  Thomson  in  1904. 10  In  his  paper,  Thomson  was  looking  to  determine  equilibrium 
configurations  of  electrons  constrained  to  the  surface  of  a  sphere  and  subject  to  Coulomb’s 
inverse-square  law.  This  problem  became  known  as  Thomson’s  problem.  To  this  day,  although 
the  solutions  to  a  handful  of  special  cases  are  known,  there  is  no  general  solution  for  n  electrons 
to  Thomson’s  problem. 

In  1942,  Laszlo  Fejes  Toth  began  studying  a  similar  problem  of  maximizing  the  minimum  pair¬ 
wise  distance  between  n  points  on  the  unit  sphere.11  This  problem  became  known  as  Fejes  Toth’s 
problem.  In  1943,  Fejes  Toth  proved  that  for  any  n  points,  there  will  always  exist  2  points  such 
that  the  distance,  d,  between  them  satisfies  the  inequality 


Furthennore,  Fejes  Toth  showed  that  this  limit  is  exact  for  n  =  3,  4,  6,  and  12;  but  like 
Thomson’s  problem,  Fejes  Toth’s  problem  has  no  known  general  solution. 

So  here  we  have  presented  2  approaches  to  the  problem  of  finding  an  “even  spacing”  of  n  points 
on  the  surface  of  the  sphere;  but  as  Whyte  (1952)  points  out,  Thomson’s  problem  and  Fejes 
Toth’s  problems  do  not  have  the  same  solution  for  general  n.  As  in  2-D  analogue  of  the  previous 
section,  we  now  state  that  the  algorithms  presented  in  this  section  are  not  exact,  or  even  optimal, 
solutions  for  producing  evenly  spaced  distributions  of  points  on  the  surface  of  the  sphere.  This  is 
a  direct  result  of  both  the  fact  that  there  is  more  than  one  formulation  of  the  problem  of  “evenly 
spacing”  n  points  on  the  surface  of  the  sphere  (we  presented  2  formulations  that  are  known  to  not 
be  equivalent),  and  that  there  was  still  no  known  general  solution  to  any  such  formulation  at  the 
time  of  this  report. 

In  fact,  only  the  vertices  of  the  5  Platonic  solids  can  be  said  to  be  truly  evenly  spaced  around  the 
surface  of  a  sphere.  A  Platonic  solid  is  any  regular,  convex  polyhedron  with  congruent  faces  of 
regular  polygons  and  the  same  number  of  faces  meeting  at  each  vertex.  Therefore,  vertices  of  the 
Platonic  solids  are  perfectly  evenly  spaced12  by  definition  (Fig.  26).  However,  there  are  only  5 
such  convex  polyhedra.13 
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Platonic  Solids 


Tetrahedron 
4  faces 


Hexahedron  (Cube) 
6  faces 


Octahedron 
8  faces 


Dodecahedron 
1 2  faces 


Icosahedron 
20  faces 


Fig.  26  The  5  Platonic  solids,  i.e.,  the  only  truly  evenly  spaced  point  distributions 
over  the  surface  of  the  sphere 


It  is  for  these  reasons  that  we  again  forego  formally  defining  the  concept  of  what  it  means  to  be 
an  evenly  spaced  distribution.  Instead,  we  offer  another  “visual  definition”  in  the  fonn  of  2 
graphs  in  Fig.  27.  The  points  in  the  graph  on  the  left  are  not  evenly  spaced,  while  those  in  the 
graph  on  the  right  are  evenly  spaced. 
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Visual  Definition  of  Even  Spacing 


Fig.  27  Visual  definition  of  what  it  means  for  points  to  be  “evenly  spaced”.  The  points  on  the  surface  of 
the  sphere  on  the  left  are  not  evenly  spaced.  The  points  on  the  surface  of  the  sphere  on  the  right 
are  evenly  spaced. 

When  evenly  spacing  points  over  the  surface  of  the  sphere,  3  of  the  most  common  techniques  are 
electrostatic  repulsion,  geodesic  subdivision,  and  spiral  methods.  Electrostatic  repulsion  is 
essentially  Thomson’s  problem;  place  n  points  anywhere  on  the  surface  of  the  sphere,  simulate 
the  repelling  forces  between  the  points,  and  look  for  a  stable  (i.e.,  converged)  configuration. 

With  geodesic  subdivision,  start  with  an  octahedron  (8  congruent  triangular  faces)  or 
icosahedron  (20  congruent  triangular  faces),  place  a  point  at  the  midpoint  of  each  edge  and 
normalize  the  coordinates  of  each  new  point  to  “push”  it  out  to  the  surface  of  the  sphere,  and 
keep  repeating  until  the  desired  number  of  points  is  attained.  Each  geodesic  subdivision 
transforms  each  triangular  face  into  4  triangular  faces,  but  as  it  involves  a  projection  to  the 
surface  of  the  sphere,  it  also  breaks  the  congruency  of  the  faces.  An  interesting  application  of 
geodesic  subdivision  is  the  design  of  domes,  buildings,  and  structures  (e.g.,  Spaceship  Earth14  at 
Epcot,  Walt  Disney  World;  the  Long  Island  Green  Dome15  in  Calverton,  NY;  the  Montreal 
Biosphere16  in  the  former  pavilion  of  the  United  States  for  the  1967  World  Fair  Expo  67  at  Parc 
Jean-Drapeau  in  Montreal,  Quebec;  and  the  Climatron17  at  the  Missouri  Botanical  Garden  in  St. 
Louis,  MO).  Geodesic  subdivision  only  works  for  configurations  of  points  where 


n  —  2  +  4k,  k  E7L+ 

(73) 

n  =  i(  4  +  5*22/c),  kET+. 

(74) 
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Configurations  with  the  number  of  points  given  by  Eq.  73  are  created  by  starting  geodesic 
subdivision  with  an  octahedron.  Configurations  with  the  number  of  points  given  by  Eq.  74  are 
created  by  starting  geodesic  subdivision  with  an  icosahedron. 

As  in  the  2-D  analogue,  spirals  are  a  natural  structure  for  distributing  points  over  the  surface  in 
question.  Just  as  electrostatic  repulsion,  geodesic  subdivision,  and  spiral  methods  are  not  the  only 
3  methods  for  distributing  points  on  the  surface  of  the  sphere,  the  following  spiral  methods  are 
neither  an  exhaustive  set  nor  do  they  claim  to  be  the  “best”  methods.  Instead,  the  spiral  methods 
presented  were  chosen  for  their  balance  of  ease  of  use  and  overall  performance. 


6.1  Rakhmanov,  Saff,  and  Zhou 

In  1994,  Rakhmanov,  Saff,  and  Zhou  were  concerned  with  what  is  commonly  referred  to  as  the 
extremal  energy  for  n  points  on  the  sphere.  To  discuss  the  extremal  energy,  we  first  need  some 
notation  and  definitions.  From  this  point  on,  all  calculations  are  for  points  on  the  unit  sphere. 
Scaling  the  points  to  spheres  of  different  radii  is  straightforward  and  left  to  the  reader. 


For  n  e  TL,  n  >  2,  let  oon  —  {xl7 ... ,  xn }  be  a  set  of  n  points  on  the  unit  sphere 

HSi  :=  {( x,y,z )  £  I3  |  x2  +  y2  +  z2  =  1}. 

Then  for  each  a  e  R,  the  a-energy  associated  with  u>n  is  defined  by 

fl!l<t<;<n  |  |  <  oc  =  0 

£(a,mn):=j  'l  a 

y2l<i<;'<n|^t  —  -fi'  |  ,  CC  A  0 

and  the  extremal  energy  for  n  points  on  is  defined  by 

£{a,n )  := 


inf  E(a,  u>n) ,  a  <  0 

(i)ncnsj 

sup  E{a,  u>n) ,  a  >  0 

a)ncfl  sf 


(75) 


(76) 


(77) 


The  problem  of  determining  the  extremal  energy  can  be  thought  of  a  generalization  of  both 
Thomson’s  and  Feje  Toth’s  problems.  In  this  notation,  Thomson’s  problem  seeks  to  detennine 
the  extremal  energy  when  a  —  —1,  i.e.,  seeks  to  detennine  £(— l,n);  Feje  Toth’s  problem  seeks 
to  determine  the  extremal  energy  when  a  —  1  [i.e.,  seeks  to  detennine  £(1,  n)]. 

Rakhmanov,  Saff,  and  Zhou  were  not  specifically  interested  in  any  one  value  of  a,  but  were 
looking  to  provide  bounds  for  the  discrete  extremal  energy  for  —2  <  a  <  2.  Furthermore,  they 
sought  a  simple  explicit  formula  for  distributing  n  points  on  the  surface  of  the  sphere  that  would 
closely  estimate  £(a,  n). 

The  reader  is  directed  to  their  1994  paper18  for  the  details  of  these  bounds,  while  the  equations  for 
what  Rakhmanov,  Saff,  and  Zhou  called  generalized  spiral  points  are  presented  in  the  following 
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equations.  The  6  and  (p  variables  have  been  switched  from  that  presented  in  Rakhmanov,  Saff,  and 
Zhou  to  confonn  to  spherical  coordinate  notation  used  throughout  this  report. 

For  any  neTL,  define  the  sequence 


hk:=l-2-^,l<k<n.  (78) 

Then  con  —  {(9k,  q)k)}k=1  c  [0,  2n]  x  [0,  n]  denotes  the  generalized  spiral  on  flS^  in  spherical 
coordinates  (p  —  1),  and  is  defined  as 

01==0n:=O,  (79) 


Qu  :  = 


6k.- 1 


mod  (2n),  2  <  k  <  n  —  1, 


(80) 


and 

(pk  :=  cos_1(hJ,  (81) 

where  C  is  a  constant  chosen  such  that  successive  points  will  all  be  approximately  the  same 
distance  apart  on  HS|.  Rakhmanov,  Saff,  and  Zhou  found  that  C  —  3.6  worked  well  for  n  <  12,000. 
Rakhmanov,  Saff,  and  Zhou’s  (RSZ)  spiral  with  1,000  points  is  shown  in  Fig.  28. 

Rakhmanov,  Saff,  and  Zhou's  Spiral 


Y-Z  View 


X-YView 


Fig.  28  One  thousand  points  of  Rakhmanov,  Saff,  and  Zhou’s  spiral  on  the  surface  of  a  sphere 

6.2  Rakhmanov,  Saff,  and  Zhou:  2  Improved  Variations 

While  an  RSZ  spiral  looks  to  cover  the  sphere  with  fairly  evenly  spaced  points,  it  can  be  easily 
seen  that  the  spacing  is  less  than  optimal  at  the  2  poles.  However,  there  is  a  very  simple 
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adjustment  that  can  be  made  to  a  standard  RSZ  spiral  to  address  the  points  in  question.  Basically, 
we  shift  the  first  and  last  points,  p0  and  pn_  1}  away  from  the  poles  and  recalculate  the  spacing  of 
the  intermediate  points.  For  the  sake  of  fast  and  efficient  implementation,  the  following 
equivalent  set  of  equations  for  an  RSZ  spiral  is  first  offered  without  proof: 


3.6 

S  ”  Vn* 


(82) 


dz  =  ^  (83) 

zk  =  1  —  k  *  dz,  0  <  k  <  n  —  1,  (84) 

rk  =  V1  -  zfe,  0  <  k  <  n  -  1,  (85) 

9k  =  k  *  s,  0  <  k  <  n  —  1,  (86) 

xk  =  rk  cos(6k),  0<k<n-l,  (87) 

yk  =  rk  sin(9k),  0  <  k  <  n  -  1.  (88) 


Then  the  RSZ  spiral  adjustment  for  addressing  the  points  at  the  poles  can  be  accomplished  by 
choosing  a  slightly  different  value  of  dz  in  Eq.  83  while  shifting  the  values  of  zk  in  Eq.  84  as 
follows: 

dz  —  -  (89) 

n  v  7 


and 

zk  =  1  —  (k  +  ^j  *  dz,  0  <  k  <  n  —  1.  (90) 

Eqs.  82  and  85-88  remain  the  same,  and  the  spiral  is  appropriately  adjusted  as  can  be  seen  in 
Fig.  29. 
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Adjusted  Rakhmanov,  Saff,  and  Zhou's  Spiral 


Y-Z  View  X-Y  View 


Fig.  29  One  thousand  points  of  adjusted  Rakhmanov,  Saff,  and  Zhou’s  spiral  on  the  surface  of  a  sphere 

Is  the  resulting  spiral  an  improvement?  That  would  depend  on  the  definition  of  what  makes  one 
distribution  of  points  better  (i.e.,  more  evenly  spaced)  than  the  other.  To  compare  the  distribution 
of  these  adjusted  RSZ  spiral  points  to  those  of  the  original  RSZ  spiral  method,  we  again  have  to 
first  establish  a  metric.  As  in  the  2-D  analogue  of  evenly  spaced  points  on  a  disc,  we  are  going  to 
calculate  Voronoi  diagrams,  but  this  time  we  will  be  working  on  the  surface  of  the  sphere  instead 
of  within  the  plane.  Our  metric  will  then  be  the  standard  deviation  of  the  surface  area  of  the  cells 
(i.e.,  of  the  solid  angles). 

Figures  30  and  3 1  show  Voronoi  diagrams  and  solid  angle  standard  deviation  for  both  the  original 
and  adjusted  RSZ  spiral  methods  for  1,000  points.  Although  both  methods  produce  distributions 
of  points  that  appear  to  be  fairly  evenly  spaced,  the  standard  deviation  of  the  solid  angle  of 
Voronoi  cells  of  the  adjusted  RSZ  spiral  method  is  slightly  more  than  41%  less  than  that  of  the 
original  RSZ  spiral  method.  This  indicates  that  the  adjusted  RSZ  spiral  method  produces  a  more 
evenly  spaced  set  of  points  as  measured  by  the  stated  metric. 
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Rakhmanov,  Saff,  and  Zhou's  Spiral 
Solid  Angle  Standard  Deviation  =  2.03871 1e-04 


Y-Z  View 


X-Y  View 


Fig.  30  Voronoi  diagram  of  Rakhmanov,  Saff,  and  Zhou’s  spiral  on  the  surface  of  a  sphere 


Adjusted  Rakhmanov,  Saff,  and  Zhou's  Spiral 
Solid  Angle  Standard  Deviation  =  1 .1 97904e-04 


Y-Z  View 


X-Y  View 


Fig.  31  Voronoi  diagram  of  adjusted  Rakhmanov,  Saff,  and  Zhou’s  spiral  on  the  surface  of  a  sphere 

Finally,  taking  a  lesson  from  Vogel’s  method  as  stated  in  Section  5.3,  we  offer  the  second 
potential  improvement  to  the  RSZ  spiral  by  changing  the  6k  increment  in  Eq.  86  from 
Rakhmanov,  Saff,  and  Zhou’s  value  of  C /yfn,  to  the  following: 


35 


8k  —  k  *  8,  0  <  k  <  n  —  1, 


(91) 


where  5  —  n{3  —  Vs)  is  the  golden  angle  as  in  Vogel’s  method.  The  other  equations  remain  the 
same  as  those  used  in  the  adjusted  RSZ  spiral  (Eqs.  82,  85,  89,  and  90)  and  the  spiral  is  re-spaced 
as  shown  in  Fig.  32. 

Adjusted  Rakhmanov,  Saff,  and  Zhou's  Spiral  with  Golden  Angle 


Y-Z  View 


X-Y  View 


Fig.  32  One  thousand  points  of  adjusted  Rakhmanov,  Saff,  and  Zhou’s  spiral  with  the  golden  angle  on 
the  surface  of  a  sphere 

Figures  3 1  and  33  show  Voronoi  diagrams  and  solid  angle  standard  deviation  for  both  the 
adjusted  RSZ  spiral  method  and  adjusted  RSZ  spiral  method  with  golden  angle  for  1,000  points. 
Although  both  methods  produce  distributions  of  points  that  appear  to  be  fairly  evenly  spaced,  the 
standard  deviation  of  the  solid  angle  of  Voronoi  cells  of  the  adjusted  RSZ  spiral  method  with 
golden  angle  is  slightly  more  than  27%  less  than  that  of  the  adjusted  RSZ  spiral  method.  This 
indicates  that  the  adjusted  RSZ  spiral  method  with  golden  angle  produces  a  more  evenly  spaced 
set  of  points  as  measured  by  the  stated  metric. 
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Adjusted  Rakhmanov,  Saff,  and  Zhou's  Spiral  with  Golden  Angle 
Solid  Angle  Standard  Deviation  =  8.70371 3e-05 


Y-Z  View  X-Y  View 


Fig.  33  Voronoi  diagram  of  Rakhmanov,  Saff,  and  Zhou’s  spiral  with  the  golden  angle  on  the  surface 
of  a  sphere 


6.3  Bauer 

The  final  method  to  be  presented  is  that  of  Bauer.19  In  2000,  Bauer  presented  a  straightforward, 
nonrecursive  spherical  spiral  for  testing  algorithms  for  stellar  attitude  determination  analyses. 
Like  Rakhmanov,  Saff,  and  Zhou’s  spiral,  Bauer’s  is  an  analytic  spiral  that  is  simple  to 
implement.  Bauer’s  spiral  has  the  form 


L  —  y/n  n, 

(92) 

2k  — 1 

zk  =  1 - , 

*  n 

1  <  k  <  n, 

(93) 

(pk  =  COS"1^), 

1  <  k  <n, 

(94) 

6k  L  (Pk> 

1  <  k  <n, 

(95) 

xk  =  sin(< pk)  cos(Ok), 

1  <  k  <  n, 

(96) 

yk  =  sin(< pk)  sin(0j, 

1  <  k  <n. 

(97) 

Bauer’s  spiral  with  1,000  points  is  shown  in  Fig.  34.  Figure  35  shows  the  Voronoi  diagram  and 
solid  angle  standard  deviation  for  the  same  set  of  points.  Here  we  can  see  that  the  standard 
deviation  of  the  solid  angle  of  Voronoi  cells  of  Bauer’s  spiral  method  is  extremely  close  to 
(about  1%  higher  than)  that  of  the  adjusted  RSZ  spiral  method  with  golden  angle.  This  indicates 
that  Bauer’s  method  and  the  adjusted  RSZ  spiral  method  with  golden  angle  produce  equivalently 
evenly  spaced  point  distributions  as  measured  by  the  stated  metric. 
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Bauer's  Spiral 


Y-Z  View  X-Y  View 


Fig.  34  One  thousand  points  of  Bauer’s  spiral  on  the  surface  of  a  sphere 


Bauer's  Spiral 

Solid  Angle  Standard  Deviation  =  8.794582e-05 


Y-Z  View 


X-Y  View 


Fig.  35  Voronoi  diagram  of  Bauer’s  spiral  on  the  surface  of  a  sphere 
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7.  Conclusions 


The  algorithms  presented  in  this  report  offer  only  a  small  glimpse  into  the  complex  and  unsolved 
problem  of  point  picking  and  distributing.  They  are  the  algorithms  that  have  proved  most  useful 
for  myself  in  a  large  variety  of  applications.  Many  of  the  algorithms  are  well  known  (e.g.,  any  of 
the  point  picking  algorithms)  but  not  always  well  documented.  In  such  cases,  it  was  my  aim  not 
only  to  document,  but  also  to  offer  some  explanation  or  proof  for  why  the  algorithm  is  correct 
and  why  it  works. 

All  but  one  of  the  point  distributing  algorithms  presented  are  formally  attributed  to  the  person  or 
people  who  discovered/derived  them.  While  I  do  not  claim  to  be  the  first  person  to  derive  what  I 
have  called  the  spiral  method,  I  was  not  able  to  find  such  a  method  and/or  equations  via  an 
extensive  literature  search  which  included  both  the  internet  and  many  scientific  and 
mathematical  journals  and  textbooks.  I  would  be  very  interested  to  know  whether  or  not  such 
equations  have  been  published  elsewhere  and  welcome  any  such  information  and/or  comments. 
Furthermore,  I  highly  suggest  that  the  interested  reader  read  through  the  sources  referenced  in 
this  report. 


39 


8.  References  and  Notes 


1.  (U)  Yager  RJ,  Patterson  C.  Modular-munition,  dual-warhead  concept  analysis.  Aberdeen 
Proving  Ground  (MD):  Army  Research  Laboratory  (US);  2013  July.  Report  No.:  ARL-TR- 
6516,  (SECRET). 

2.  (U)  Strohm  L,  Arthur  M.  Modular  lethality:  examining  changes  in  module  placement  for 
multi-target  defeat.  Aberdeen  Proving  Ground  (MD):  Army  Research  Laboratory  (US);  2015 
May.  Report  No.:  ARL-TR-7246,  (SECRET). 

3.  Strohm  L.  Modular  lethality:  methodology  for  resource  allocation.  Aberdeen  Proving 
Ground  (MD):  Army  Research  Laboratory  (US);  in  press. 

4.  Muller  ME.  Some  continuous  Monte  Carlo  methods  for  the  Dirichlet  problem.  Ann  Math 
Stat.  1956;27(3):569-589. 

5.  Muller  ME.  A  note  on  a  method  for  generating  points  unifonnly  on  n-dimensional  spheres. 
Communications  of  the  ACM.  1959;2(4):  19-20.  Note:  There  is  an  error  in  the  sign  of  the 
exponent  in  the  equation  for  step  2  of  Muller’s  Method  in  this  reference. 

6.  Gautschi  W.  The  spiral  of  Theodorus,  numerical  analysis,  and  special  functions.  Journal  of 
Computational  and  Applied  Mathematics.  2010;235(4):  1042-1052. 

7.  Vogel  H.  A  better  way  to  construct  the  sunflower  head.  Mathematical  Biosciences. 
1978;44(3-4):  178-189. 

8.  Vogel  referred  to  the  spiral  as  a  cyclotron  spiral,  but  the  author  was  unable  to  find  any  other 
resources  that  refer  to  the  spiral  in  question  by  this  name. 

9.  Vogel’s  method  is  actually  a  single  spiral  method. 

10.  Thomson  JJ.  On  the  structure  of  the  atom:  an  investigation  of  the  stability  and  periods  of 
oscillation  of  a  number  of  corpuscles  arranged  at  equal  intervals  around  the  circumference  of 
a  circle;  with  application  of  the  results  to  the  theory  of  atomic  structure.  Philosophical 
Magazine.  1904;7(39):237-265.  Series  6. 

1 1 .  Whyte  LL.  Unique  arrangements  of  points  on  a  sphere.  The  American  Mathematical 
Monthly.  1952;59(9):606-61 1. 

12.  All  points  are  equidistant  from  their  nearest  neighbors. 

13.  Proofs  can  be  found  in  a  variety  of  reference  books.  For  one  such  proof,  see  section  4.1.2  of 
O-Rourke  J.  Computational  Geometry  in  C.  Cambridge  (United  Kingdom):  Cambridge 
University  Press;  1993. 


40 


14.  Spaceship  Earth.  Walt  Disney  World;  [accessed  2014  Nov  7] 
https://disneyworld.  disney .  go  .com/ attractions/ epcot/ spaceship-earth/. 

15.  Long  Island  Green  Dome,  [accessed  2014  Jul  7] 
https://www.facebook.com/ligreendome/map?activecategory=Photos&session_id=1334292191. 

16.  Environment  Canada.  Biosphere,  environment  museum,  [accessed  2014  Nov  7] 
http://www.ec.gc.ca/biosphere/. 

17.  Missouri  Botanical  Garden.  Climatron.  [accessed  2014  Nov  7] 
http://www.missouribotanicalgarden.org/gardens-gardening/our-garden/gardens- 
conservatories/conservatories/climatron.  aspx . 

18.  Rakhmanov  EA,  Saff  EB,  Zhou  YM.  Minimal  discrete  energy  on  the  sphere.  Mathematical 
Research  Letters.  1994;l(6):647-662. 

19.  Bauer  R.  Distribution  of  points  on  a  sphere  with  application  to  star  catalogs.  Journal  of 
Guidance,  Control,  and  Dynamics.  2000;23(  1):  1 30—137. 


41 


INTENTIONALLY  LEFT  BLANK 


42 


Appendix.  Limit  Evaluation  from  Section  5.1 


This  appendix  appears  in  its  original  form,  without  editorial  change. 
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Claim:  lim  lim  2k  +  1  —  2 4k  Vk  +  1  cos  f^k+1  =  -3— 

n-> oo  k->n  \  b  J  4 

Proof:  Let  /(x)  =  2x  +  1  —  2 Vx  Vx  +  1  cos  ^  Then 
lim  /(x)  =  lim  1  +  2x  —  2 Vx(x  +  1)  cos  /^*+1 

x— >oo  x— >oo  V  b  / 

=  lim  i  +  2x  -  zV^Tlj  cos  ( f=±§)) 

=  lm  1  +  2x  -  2,/*0  +  1)  cos 
=  lim  1  +  2x  -  2^  +  1)  cos  (i  (^g))  • 


Let  z  —  — .  Then 

X 


lim  /(x)  -  lim/ 

x->oo  z->0 


=  lim  1  +  -  —  2 

z->0  z 


Let  g(z)  =  J“  (“  +  ar|d  calculate  the  Laurent  series  about  z  —  0. 


g(z)  =  J“(z  +  -0  =  Sn=-ooanZr 


where 


an  =  —  cf> 

9  m  ■'1 


S(z) 


27TI  y  Z 


dz. 


n+l  ' 


and  the  path  of  integration,  y,  is  counterclockwise  around  any  closed,  rectifiable  path  containing 
no  self-intersections,  and  enclosing,  but  not  containing,  the  point  z  =  0. 
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Note  that 


an  =  —$>  °g-dz  =  ±$ 

n  2m  tY  zn+1  2m  Ty 


m+i) 


dz  —  —S  -j(l  +  z)dz 
2m  J Y  y 


=  — <f>  -777 V 1  +  z  dz  =  —6  A(z)dz 

2ni  rY  zn+2  2ni  rY  y  J 

where  A(z)  —  -777VI  +  z,  and  that  A(z)  is  analytic  everywhere  for  n  <  —2. 
Therefore,  by  the  Cauchy  integral  theorem* 
an  —  —  <f>  A(z)dz  =  0  forn  <  —2. 

n  2m  JY 

So  we  need  only  to  consider  n  =  —  1, 0, 1, 2, ... 
n  =  —1: 


a_i  =  —  <f>  -V 1  +  zdz 

1  2m  JY  z 

z  —  0  is  a  pole  of  order  1 ,  therefore 

a-i  =  h,i  =  Res» (;Vl+7)  =  to (z gVl+i))  =  1. 

n  —  0: 

a0  —  —  <5  — 7  Vl  +  zdz 

u  2m  JY  Z 2 

z  =  0  is  a  pole  of  order  2,  therefore 

“»  =  hi  ^'/T+Idz  =  Res0  (iVTTz)  =  to^(z2  (pVTTl)) 
=  lim-(l  +  z)-1/ 2  =  i 

z->0  2  2 


n  =  1: 


a-.  =  — <f>  -7 Vl  +  zdz 

1  27Tt  Ty  z3 


0  is  a  pole  of  order  3,  therefore 

a^  =  hi  7VTTI<iz  =  Res«  (7VITz) 

=  -lim-^fz3  (-^Vl  +  z))  =  -lim  — -(1  +  z)~3/ 2  = 

2  z^o  dz2  V  Vz3  7/  2  z^o  4  v  ' 


Also  know  as  the  Cauchy-Goursat  theorem:  Let  D  Q  C,D  open  and  simply  connected,  let  f:D  -*  C  be  a  holomorphic 
(analytic)  function,  and  let  y  be  a  rectifiable  path  in  D  whos  start  point  is  equal  to  its  end  point.  Then  $  f(z)cLz  =  0. 
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Therefore 


s&  =  J;G+1)  =  a^n  =  r +  r - 1  ■ + 


Vz 
.1+Vz+l^ 


fi. 


Now  let  h(z )  =  cos  ^  (1+^+1)^  and  since  h(z)  is  differentiable  at  z  =  0,  the  Maclaurin 
of  /i(z)  is  given  by 

h(z)  =  h(0)  +  h'(0)z  +  0(z2). 

Therefore, 

h'(z)  =  ( - — - 7 - 1  ,  .■)  sin  f- 

V2dVJ+i(i+v5+i)  2ftVi(i+v^+i);  yA: 

So  let  h'(z)  =  hx(z)  —  h2(z)  where 

/^(z)  =  ( - — - -)  sin  ^ 

and  h*  ^  =  (iw^vSf)) sln  (;  (irlff))- 

h2  (z)  is  not  defined  at  z  =  0,  but  taking  the  limit  as  z  -»  0  yields 

lim  h2 (z)  —  lim  ( — 1  _ - )  sin  f  - ( — ^==] ) 

/v  y  \2i>Vz(l+Vz+l) /  \  ft  Vl+Vz+1/  I 


series 


sm(i>(  i+3+i)) 

—  lim — A — t — 

z->0  2fcVz(l+Vz+l) 

lim  sin(^-( — 

z^o  y£>  v  l+vz+i/y  o 

lim  2i>Vz(l+Vz+l)  o’ 

and  applying  L’ Hospital’s  rule  yields 


r  .  dz(sin(fo(i+^rr) )) 

lim  n2  (z)  -  lim  -3-7 — —  ( 

z->0  z->0  — ^21>Vz(l+Vz+l)) 


cos 


..  IbU+Vz+l//  1 

_ lim _ _ _ _ _ _ — 

z->0  fe(l+V’z+l)(l+2z+Vz+l)  8b2’ 


So  the  singularity  of  h2  (z)  at  z  =  0  can  be  removed,  and  h2  (z)  becomes 
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hz  (z)  (Um) sin  (r  Griff))  >  °therwise 

Therefore 

h'(z)  =  hi(z)  -  h2(z) 

^h'(0)  =  h1(0)-/i2(0)  =  0-i^ 

=>  ^0)  =  C0S(Kfflff))  =  ^(°)  +  h'(0)z  +  O(z2)  =  1  -  Jl  +  °(z2)- 

Finally,  ff(z)  =  ,];(;  +  l)  =  ;  +  +  0(z2) 

and  /t(z)  =  coS^(TT^=))  =  l-j|T  +  0(z2) 
imply  that  lim  /(x)  =  lim/  (-) 

x->oo  z->  0  Vz/ 

=  >™0 1 + ; "  2J;(;+1)  cos (f (iris)) 

—  lim  1  +  -  —  2  g(z)  h(z) 

z->0  z 

=  !™„1  +  ;-2  (: + 1 _  I +  0(22:i)  (x  _  +  °^) 

=  i1?o1+;-2  (;+;-sp+0(z)) 

=  lim  1  +  -  -  -  -  1  +  +  0  (z)  = 

z->0  z  z  4£>2  v  y  4£>2 
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